% This file generates the IRF presented in the sixth quadrant (2,3) in 
% Figure 10 in the paper. It needs to be ran from 
% 'FIG_10_Robustness_Inflation.m'.

global P2 lambda P H 
%% Options for the specific results to generate.

graph           =   0; 

csv_file        =   1;          % = 1 if you want to import from CSV.
                                % = 0 if you want to import from XLSX.
                                
use_data        =   1; 
                                % = 1  if inflation using Michigan Expectations
                                % Median. BASELINE    
                                % = 2  if inflation using Michigan Expectations
                                % Mean.                     
                                % = 3  if inflation using SPF CPI Median.
                                % = 4  if inflation using SPF GDP Deflator.
                                % ahead
                                % = 5 if inflation using the Cleveland Expectations.
                                
%% Importing and Preparing the data.

if csv_file == 1
    data            =   csvread('Data_File_Baseline.csv');        
else
    data            =   xlsread('Data_File.xlsx','Baseline');
end
% Import the specific variables. 
dz          =   data(:,3);
dy          =   data(:,4);
zlb         =   data(:,5);
rec         =   data(:,6);
zlev        =   data(:,7);
yy          =   data(:,8);
time        =   data(:,9);
    
if use_data == 1        % Michigan Survey of Consumer Expectations Median.
    ylev_data   =   csvread('Data_File_Inflation.csv');
    ylev        =   ylev_data(:,6);   
elseif use_data == 2    % Michigan Survey of Consumer Expectations Mean.
    ylev_data   =   csvread('Data_File_Inflation.csv');
    ylev        =   ylev_data(:,7);   
elseif use_data == 3    % SPF forecast of CPI inflation.
    ylev_data   =   csvread('Data_File_Inflation.csv');
    ylev        =   ylev_data(:,8);   
elseif use_data == 4    % SPF forecast of GDP deflator inflation.
    ylev_data   =   csvread('Data_File_Inflation.csv');
    ylev        =   ylev_data(:,10);    
elseif use_data == 5    % Cleveland's Fed expected inflation.
    ylev_data   =   csvread('Data_File_Inflation.csv');
    ylev        =   ylev_data(:,18); 
end

% Defining the sample:
st          =   148;                        % = 1 for full sample 
                                            % = 148 for 1984Q1 
                                            % = 52 for 1960Q1.  
                                            
%% Main loop to generate IRFs.

Ps              =   [2; 3; 4];

for ix = 1:length(Ps)
    
    P                       =   Ps(ix);
    
    % LHS:
    ylev2       =   (1/4)*ylev(st:end);
    % Productivity:
    zlev2       =   zlev(st:end);
    % Dummy for ZLB:
    zlb2        =   zlb(st:end);
    % Inflation: 
    yy2         =   yy(st:end);
    % Creating W matrix: 
    ff          =   [ylev2 zlev2 yy2];
    w           =   lagmatrix(ff,1:P);  % Remove period t observations.

    % Creating lagged variables interacted with period-t ZLB dummy:
    [TY,YT]     =   size(w);
    zlb3        =   zeros(TY,YT);
    for jx = 1:YT
        zlb3(:,jx)  =   zlb2;
    end
    www         =   zlb3.*w; 
    zlb4        =   1-zlb3;
    wwww    	=   zlb4.*w;

   % LHS variable:   
    y                       =   ylev2(1+P:end,:); 

    % Local projection.
            % Standard case.  
    w                       =   [zlb2(1+P:end) zlb2(1+P:end).*zlev2(1+P:end),...
                                 wwww(1+P:end,:), www(1+P:end,:)];
    x                       =   (1-zlb2(1+P:end)).*zlev2(1+P:end,:); % Shock variable.   
    lp                      =   locproj_var(y,x,w,h1,H,'reg'); 
    % Extract the IRF from local projection: 
    irf                     =   lp.IR(2:H+1); 
    % Construct standard errors:
    u                       =   lp.Y - lp.X*lp.theta;
    [nwse_lp,vc_lp]         =   NeweyWest(u,lp.X,0,0);
        % ZLB case.
    w                       =   [zlb2(1+P:end) (1-zlb2(1+P:end)).*zlev2(1+P:end),...
                                 wwww(1+P:end,:), www(1+P:end,:)];
    x                       =   zlb2(1+P:end).*zlev2(1+P:end,:);  
    lp_zlb                  =   locproj_var(y,x,w,h1,H,'reg'); 
    % Extract the IRF from local projection: 
    irf_zlb                 =   lp_zlb.IR(2:H+1); 
    % Construct standard errors:
    u                       =   lp_zlb.Y - lp_zlb.X*lp_zlb.theta;
    [nwse_lp_zlb,vc_lp_zlb] =   NeweyWest(u,lp_zlb.X,0,0);

    % Smooth local projection.
            % Standard case.
    w                       =   [zlb2(1+P:end) zlb2(1+P:end).*zlev2(1+P:end),...
                                 wwww(1+P:end,:), www(1+P:end,:)];
    x                       =   (1-zlb2(1+P:end)).*zlev2(1+P:end,:); 
    slp                     =   locproj_var(y,x,w,h1,H,'smooth',r,lambda); 
    % Extract the IRF from local projection: 
    irf_slp                 =   slp.IR(2:H+1);
    P2                      =   slp.P2;
    % Construct standard errors:
    u                       =   slp.Y - slp.X*slp.theta;
    [nwse_slp,vc_slp]       =   NeweyWest_slp(u,slp.X,0,0);
    V2                      =   slp.B*vc_slp(1:slp.K,1:slp.K)*slp.B';
    se_slp                  =   sqrt(diag(V2));

            % ZLB case.
    w                       =   [zlb2(1+P:end) (1-zlb2(1+P:end)).*zlev2(1+P:end),...
                                 wwww(1+P:end,:), www(1+P:end,:)];
    x                       =   zlb2(1+P:end).*zlev2(1+P:end,:); 
    slp_zlb                 =   locproj_var(y,x,w,h1,H,'smooth',r,lambda); 
    % Extract the IRF from local projection: 
    irf_slp_zlb             =   slp_zlb.IR(2:H+1);
    % Construct standard errors:
    u                       =   slp_zlb.Y - slp_zlb.X*slp_zlb.theta;
    [nwse_slp_zlb,vc_slp_zlb] = NeweyWest_slp(u,slp_zlb.X,0,0);
    V22                     =   slp_zlb.B*vc_slp_zlb(1:slp_zlb.K,1:slp_zlb.K)*slp_zlb.B';
    se_slp_zlb              =   sqrt(diag(V22));

   
     %% Saving results to Excel.
   
    p = 0.9;
    bb = tinv(p,TY-P);
    
    if P == 2
        irf_lp_P_2                  =   irf;
        irf_slp_P_2                 =   irf_slp;
        irf_lp_P_2_zlb              =   irf_zlb;        
        irf_slp_P_2_zlb             =   irf_slp_zlb;
        irf_slp_P_2_LB              =   (irf_slp-bb*se_slp(1:H));
        irf_slp_P_2_UB              =   (irf_slp+bb*se_slp(1:H));
        irf_lp_P_2_LB               =   (irf-bb*nwse_lp(1:H));
        irf_lp_P_2_UB               =   (irf+bb*nwse_lp(1:H));
    elseif P == 3
        irf_lp_P_3                  =   irf;        
        irf_slp_P_3                 =   irf_slp;
        irf_lp_P_3_zlb              =   irf_zlb;        
        irf_slp_P_3_zlb             =   irf_slp_zlb;
    elseif P == 4
        irf_lp_P_4                  =   irf;
        irf_slp_P_4                 =   irf_slp;
        irf_lp_P_4_zlb              =   irf_zlb;        
        irf_slp_P_4_zlb             =   irf_slp_zlb;
    end
end
